#include <iostream>
#include <fstream>
#include <string>
#include <sstream>

#include "Meshes/StructuredMesh.hpp"
#include "Meshes/Uniform.hpp"
#include "Utils/inout.hpp"

#include "ParTrack/ParticleTracking.hpp"
#include "ParTrack/Linear.hpp"
#include "ParTrack/RungeKutta4.hpp"

using namespace std;
using namespace Tuna;

typedef TunaArray<double, 2>::tiny floatTinyArray_t;
typedef TunaArray<double, 2>::huge ScalarField2D;

void writeTofile_DX(floatTinyArray_t *, int, string, int, int);

int main () 
{
    StructuredMesh< Uniform<double, 2> > mesh(1, 257, 1, 257);
    floatTinyArray_t deltas = mesh.getDeltas();
    double dt = 0.0001;
    int Npoints_x = 25, Npoints_y = 25, Nsteps = 10, frec = 1, Ncicles = 50;
    
    floatTinyArray_t plane[Npoints_x * 2 * Npoints_y];

    double left = 0.1;
    double right = 0.9;
    double top = 0.49;
    double bot = 0.1;

    double factor_x = (right - left) / (Npoints_x - 1);
    double factor_y = (top - bot) / (Npoints_y - 1);

    cout << "\n ----- Construction of two planes ----- \n\n";
    
    for(int j = 0; j < Npoints_y; ++j) 
	for(int i = 0; i < Npoints_x; ++i) {
	    plane[i + j * Npoints_x](0) = left + i * factor_x;
	    plane[i + j * Npoints_x](1) = bot  + j * factor_y;
	}
    bot = 0.51;
    for(int j = Npoints_y; j < 2 * Npoints_y; ++j) 
	for(int i = 0; i < Npoints_x; ++i) {
	    plane[i + j * Npoints_x](0) = left + i * factor_x;
	    plane[i + j * Npoints_x](1) = bot  + (j - Npoints_y) * factor_y;
	}

    writeTofile_DX(plane, 0, "./PlaneMixing/plane.", Npoints_x, Npoints_y);

    cout << "\n ----- Tracking: begins ... \n\n";

// Velocity field at time t
    ScalarField2D u(mesh.getExtentNodes());       
    ScalarField2D v(mesh.getExtentNodes());

// Velocity field at time t + dt
    ScalarField2D un(mesh.getExtentNodes());       
    ScalarField2D vn(mesh.getExtentNodes());

    ParticleTracking< RungeKutta4< Linear<double, 2> > > particle;
    
    int count = 1;
    for(int cicles = 1; cicles <= Ncicles; ++cicles) {
	cout << " ----> Cicle " << cicles << endl;
	for(int step = 1; step <= Nsteps; ++step) {
	    cout << " step = " << step 
		 << " count = " << count << endl;

// Read the velocity at time t from a file
	    InOut::readFromFile_1(u, v, step, "./Velocity/velc.");
// Read the velocity at time t + dt from a file
	    InOut::readFromFile_1(un, vn, step + 1, "./Velocity/velc.");

// Track each particle on the plane
	    for(int j = 0; j < Npoints_y; ++j) 
		for(int i = 0; i < Npoints_x; ++i)
		    plane[i + j * Npoints_x] = 
			particle.track( plane[i + j * Npoints_x], 
					deltas, u, v, un, vn, dt );

	    for(int j = Npoints_y; j < 2*Npoints_y; ++j) 
		for(int i = 0; i < Npoints_x; ++i) {
		    plane[i + j * Npoints_x] = 
			particle.track( plane[i + j * Npoints_x], 
					deltas, u, v, un, vn, dt );
		}
	    
	    if( !(count % frec) )
		writeTofile_DX(plane, count, "./PlaneMixing/plane.", Npoints_x, Npoints_y);
	    count++;
	}
    }
    
    cout << "\n ----- Tracking: ends ----- \n\n";
    return 0;
}

void writeTofile_DX(floatTinyArray_t *plane, int n, string filename,
		    int Npoints_x, int Npoints_y)
{
    ostringstream extnumber;
    extnumber.width(5);
    extnumber.fill('0');
    extnumber << n;
    string snumber = extnumber.str();
    filename += snumber;
    ofstream file (filename.c_str());  

    cout << "writing " << filename << endl;
    
    file << " object 1 class array type float rank 1 shape 2 items "
	 << Npoints_x * 2 * Npoints_y << " data follows\n";
    for(int j = 0; j < 2 * Npoints_y; ++j) 
	for(int i = 0; i < Npoints_x; ++i)
	    file << plane[i + j * Npoints_x](0) << "  " 
		 << plane[i + j * Npoints_x](1) << "  ";
    file << endl;

    file << "object 2 class array type int rank 1 shape 2 items "
	 << 2 << " data follows\n";
    for(int i = 0; i <  2; ++i)
	file << i << "   " << i + 1 << "   ";
    file << endl;
    file << "attribute \"element type\" string \"lines\"\n"
	 << "attribute \"ref\" string \"positions\"\n";

    file << "object 3 class array type float rank 0 items "
	 << Npoints_x * 2 * Npoints_y << " data follows\n";

    for(int i = 0; i < Npoints_x; ++i)
	for(int j = 0; j < Npoints_y; ++j) 
	    file << 1.00 << "  ";
    for(int i = 0; i < Npoints_x; ++i)
	for(int j = Npoints_y; j < 2*Npoints_y; ++j) 
	    file << 1.01 << "  ";
    file << endl;
 
    file << "attribute \"dep\" string \"positions\"\n";
    file << "object \"irregular positions regular connections\" class field\n"
	 << "component \"positions\" value 1\n"
	 << "component \"connections\" value 2\n" 
	 << "component \"data\" value 3\n" 
	 << "end\n" ;
    file.close();
}


